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Abstract 

We consider the optimal spacing between marine reserves for maximising the vi- 
abihty of a species occupying a reserve network. The closer the networks are placed 
together, the higher the probability of colonisation of an empty reserve by an occu- 
pied reserve, thus increasing population viability. However, the closer the networks 
are placed together, the higher the probability that a catastrophe will cause extinction 
of the species in both reserves, thus decreasing population viability. Using a simple 
discrete-time Markov chain model for the presence or absence of the species in each 
reserve we determine the distance between the two reserves which provides the optimal 
trade-off between these processes, resulting in maximum viability of the species. 
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1 Introduction 

The design of reserve networks for marine conservation is a contentious issue from 
both theoretical and practical points of view. One of the difficult theoretical issues 
in marine reserve system design concerns the optimal connectedness of the reserves 
within the network. The closer two reserves are placed together, the more likely it is 
that a population occupying one reserve will colonise the other, increasing population 
viability. However, most theories fail to consider the possibility of catastrophic events 
(Possingham et al. (2000), Shafer (2001), Sala et al. (2002) and Allison et al. (2003)); 
the closer two reserves are placed together the more likely they are to be struck by the 
same catastrophic event, thus decreasing viability. Consequently, there exists a natural 
tension between these two processes. Whilst this trade-off has been identified by 
many authors, a detailed theoretical investigation of its influence and the subsequent 
analytic determination of an optimal spacing, has, to the best of our knowledge, not 
been undertaken. 
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In this paper we explore the question of how closely two marine reserves should 
be placed in order to maximise the viability of a species in a reserve network under 
the competing processes of colonisation and catastrophe. We use a simple discrete- 
time Markov chain metapopulation model to explore the effect of reserve spacing on 
the viability of the metapopulation. The mathematical fomulation of the model fol- 
lows closely the work of Day and Possingham (1995). We use the quasi-stationary 
distribution (the distribution of the metapopulation conditioned upon the population 
being extant) to determine an optimal inter-reserve spacing. The analysis relies on 
finding eigenvalues of the transition matrix of the Markov chain. The rate at which 
the metapopulation decays to extinction from quasi-equilibrium is given by the second 
eigenvalue of the probability transition matrix (Darroch and Senata (1965)). Thus, 
having found the second eigenvalue of the Markov chain, we can investigate the in- 
fluence of intcr-rcscrvc distance on metapopulation viability, and by maximising this 
eigenvalue determine the optimal inter-reserve distance. 

Having constructed a baseline model for the viability of a marine metapopulation, 
we will then examine how the inclusion of two further processes impacts optimal re- 
serve spacing. The inclusion of an external recruitment process allows us to model a 
more expansive population which is not just confined to the reserves under considera- 
tion (Sweaerer et al. (2002) and Cowen et al. (2000)). The analysis of this new model 
is similar to that for the baseline model, however we must now use the stationary dis- 
tribution to investigate the effect of external recruitment since extinction (absorption) 
is not possible. We determine the optimal spacing by minimising the probability of 
having both reserves empty in equilibrium. We extend the model further by including 
the possibility of local reserve extinction, in which the population can become extinct 
in the absence of catastrophes. These two additional processes allow us to create a 
more general and realistic model of a marine metapopulation and to develop a more 
comprehensive idea of optimal inter-reserve spacing. This is a significant extension 
over the recent work by McCarthy et al. (2005) and a related uncertainty analysis of 
that model (Halpern et al. (2006)). 

In the baseline case, and the two extensions to this model, we find reserve spac- 
ings that maximise the viability of the metapopulation. Our analytical results are 
an advance on existing theory, which is largely driven by results from Monte Carlo 
simulation models and often ignores the effect of catastrophes. 

2 The Baseline Model 

Our baseline model consists of two processes: colonisation, and extinction caused by 
catastrophic events. This model is constructed to find the optimal spacing between 
two marine reserves under these competing processes. A catastrophe in this case is an 
extinction event which affects one or both populations. This occurs as either a direct 
reduction in population size or a process which degrades the habitat to such an extent 
that no organism will survive after a given time period, e.g. due to coral bleaching 
(Alhson et al. (2003)). 

For each reserve we model the presence or absence of the marine species within 
that reserve (Akcakaya and Ginzburg (1991)). We also assume that both patches are 
identical, thus reducing the state space to three states since we only need to keep track 
of whether zero (total extinction), one or both patches are occupied. 
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Figure 1: Schematic of an example scenario where an oil spill has occurred along a coast 
line with two defined marine reserves which are placed at some distance apart. 

We have used basic probability theory to construct a discrete-time Markov chain 
for each of the processes - colonisation and catastrophe - considered in this model. To 
construct our probability transition matrix we need to determine the probabilities for 
each possible transition in our state space. Our model examines two marine reserves 
which lie at some point along a coastline. These two reserves are separated by a 
distance d and are susceptible to being struck by some catastrophic event, such as an 
oil spill, which has a probability r of occurrence in each time step. We assume that 
the probability that the catastrophe affects the second reserve decays exponentially 
with the distance d between the reserves, and additionally depends upon the mean 
catastrophe size /x. Based upon these assumptions we construct an extinction matrix 
E. An empty reserve may be colonised by an occupied reserve with a probability that 
decays exponentially with the distance d between the reserves, and also depends upon 
a fitness parameter a, which may be estimated for a particular marine system and 
accounts for the landscape between a particular reserve pair. The fitness parameter a 
may also be interpreted as the inverse of the average distance an individual disperses 
over the time period. Based upon these assumptions we construct a colonisation 
matrix C. Having constructed these matrices we can calculate the transition matrix 
A for the complete process. We assume that extinction occurs before colonisation, 
and thus A = E x C (Gilpin (1992), Day and Possingham (1995)). This order of 
matrix multiplication is usually chosen since it is a "worst-case" scenario, in that 
it will overestimate the probability of extinction. Additionally, in the present study 
we are only considering the second eigenvalue of the matrix A, and thus the order 
of multiplication does not affect the results. To calculate the probability of survival 
at some future time we may recursively evaluate Pt+i = PtA, where Pt is a vector 
containing the probabilities of being in each state at time step t. 

2.1 Colonisation Process 

We now need to establish the colonisation probability matrix. In this model we have 
that an unoccupied reserve can be colonised by an occupied reserve in each time step. 
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Initially wc construct a state space diagram (Figure (2)), which shows the four possible 
transitions. If both reserves are empty, then there are no individuals in the reserve 
network to colonise the reserves, and thus the zero state is clearly absorbing. 



Figure 2: State space with transitions between the numbers of occupied patches for the 
colonisation transition paths in our two patch metapopulation. 

The process of migration from one patch to another is dependant upon the fitness 
a of a population within the landscape and the distance d that the two reserves are 
placed apart from each other. This form of colonisation probability is ubiquitous in 
the literature and has been empirically observed for many metapopulations (Gilpin 
(1992), Hanski (1994a, 1994b) and Cowen et al. (2000)). Hence, the colonisation 
matrix is given by 



In Figure (3) we can see that the chance one patch is colonised from another decays 
with distance and the relationship is dependent upon the fitness parameter a. In a 
marine environment, the process of colonisation takes place continuously through time, 
implying that the dispersal of marine organisms in all stages of their life cycle, whether 
that is dispersal of larvae or the translocation of adults, plays an important role in 
metapopulation viability (Cowen et al. (2000)). Furthermore, organisms who are able 
to translocate as a stress response, related to some catastrophic event affecting their 
environment, will also migrate to an unaffected reserve in the manner proposed by C. 

2.2 Extinction Process 

We now derive the extinction probability matrix. Once again we construct a state 
space diagram (Figure (4)), which shows the six possible transitions. Since we are 
only considering the process of extinction, an increase in the number of occupied 
reserves is impossible. 

A catastrophe that affects at least one reserve, arrives with probability r within 
each time step. When a catastrophe occurs, we assume that it affects both reserves 
with probability exp(— d//^), and therefore affects only one reserve with probability 
1 — cxp(— (i//u), where /j, is the mean catastrophe size. This formulation assumes that 
the probability that a catastrophe affects both reserves decays exponentially with 
increasing inter-reserve distance. This is a reasonable assumption, in particular for 
catastrophes with a radius that is exponentially distributed in size, with an average 
radius of /i. 

Wc shall now derive the functional forms for each possible transition. Firstly, if 
neither reserve is occupied, then the system will remain unoccupied since we are only 
considering catastrophic events. Thus, the entry £^(0, 0) of our transition probability 





(1) 
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Figure 3: Graph of the probabihty of colonisation exp{—ad) with respect to distance d, for 
a — 0.01 (sohd hne), 0.1 (dots) and 0.2 (dashed hne). 




Figure 4: State space with transitions between the number of occupied patches for the 
catastrophe transition paths in our two patch metapopulation. 



matrix E will be 1, and £'(0, 1) and £"(0, 2) will both be 0. When both reserves are 
occupied, we have that the system will only remain fully occupied if a catastrophe does 
not occur within the time step, thus E{2, 2) = 1 — r. If a catastrophe does occur (with 
probability r), it removes only one reserve with probability 1 — exp(— d///) and both 
reserves with probability exp(— Thus, E{2, 1) and £(2,0) are r(l — exp(— d///)) 
and r exp(—d/ n), respectively. Finally, when only one reserve is occupied, we have 
that the system will remain with only one reserve occupied if either a catastrophe 
does not occur (with probability 1 — r), or a catastrophe occurs but only affects 
the unoccupied reserve - this occurs with probability (1/2)(1 — exp(— d/^)), since 
a catastrophe affects only one reserve with probability (1 — exp(—d/fi)), and with 
probability 1/2 this was the one occupied (since we are assuming that each patch is 
identical and that the catastrophe is equally likely to affect either patch initially). 
Thus, £(1,1) = (1 — r) + (r/2)(l — exp(-d/;L()). The state of the system will change to 
no reserves occupied if the catastrophe affects both reserves, or affects only one reserve, 
with this reserve being the one that was initially occupied. Therefore, £(1,0) = 
r[(l/2)(l — exp(— d/yu)) + exp(— d///)], and since we are only considering catastrophic 
events, £(1, 2) = 0. 
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This generates our transition matrix for the catastrophe/extinction process, 
/I \ 



E = 



i(l + e^ 



V 



re 



l-i(l + e- 
r ^1 — e"*^^ 





1-r 



(2) 



Figure (5) shows the probabihty r exp(—d/ fi) of a catastrophe occurring in a time 
step and it affecting both reserves, for different values of n and r = 0.5. 




Figure 5: The probability of the extinction process hitting both reserves in a network within 
a year with respect to its mean size; 11 — 20 (sohd hne) , 30 (dots) and 40 (dashed hne) , all 
with respect to the occurrence of a catastrophe in one year r — 0.5. 



3 Baseline Model Results 

As mentioned previously, to analyse the baseline model we use the second eigenvalue 
of the probability transition matrix. This eigenvalue provides the rate at which the 
metapopulation decays to extinction from quasi-stationarity. The graph below (Figure 
(6)) plots the second eigenvalue of^ = £^xCasa function of inter-reserve spacing 
d. The optimal inter-reserve distance is given by the value of d for which the second 
eigenvalue is largest. 

We can see from Figure (6) that when the mean catastrophe size jj, is smaller than, 
or comparable to, the mean dispersal distance 1/a, there exists a unique optimal inter- 
reserve spacing, signified by the hump in the curves corresponding to n = 5 and 10. 
It can also be seen that the optimal inter-reserve distance increases with increasing 
catastrophe size. As /x increases, eventually the trade-off between between colonisa- 
tion and catastrophe disappears, and the optimal inter-reserve distance is found to be 
much larger. This emerges because as the distance increases (so as to attempt to avoid 
catastrophes affecting both reserves) , the probability of colonisation becomes smaller. 
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Figure 6: Second eigenvalue of A; r — 0.5, a — 0.1 and fj, — 5 (solid line), 10 (dots), and 20 
(dashed line), as three different values for the mean size of disaster. 



Eventually colonisation becomes essentially impossible, and as a consequence the opti- 
mal inter-reserve spacing is almost entirely determined by attempting to minimise the 
probability of a catastrophe affecting both reserves. In this situation, it is still optimal 
to place the reserves at the minimum distance for which the curve is no longer increas- 
ing (since catastrophes will usually be infrequent compared to colonisation events). 
This is to allow for the best chance of recolonisation. 

Given particular parameter estimates, oTir model and method of analysis allows 
the calculation of the optimal inter-reserve distance. This generality is very helpful in 
attempting to determine the optimal placement of reserves within a complex marine 
park such as the Great Barrier Reef Marine Park, off Australia's east coast. Rather 
than having to use Monte Carlo simulations to estimate the probability of survival at 
a given distance in a two patch model such as this, we can calculate directly the effect 
of distance and plan more effectively for disaster. 

4 Model with External Recruitment 

In marine reserves, very few species have such a restricted occurrence that they only 
exist in two places. Thus, a logical extension of the model is to allow for external 
recruitment into the two patch metapopulation. This local patch effect, which has 
been coined the rescue effect by Brown and Kodric-Brown (1977), proposes that the 
migration of individuals into a struggling population may allow for the reduction in 
local extinction. The state space diagram in Figure (7), shows how the inclusion of 
external recruitment allows a locally extinct population (zero patches occupied), to 
now be recolonised (Jones et al. (1999) and Strathmann et al. (2002)). As such, 
the inclusion of external recruitment into our model creates a system without an 
absorbing state. However, we will refer to the state with both reserves unoccupied as 
being extinct. 

To construct the new overall transition matrix A we first need to construct an ex- 
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Figure 7: State space with transitions between the number of occupied patches for the 
recruitment transitions paths in our two patch metapopulation. 



ternal recruitment matrix R, that tells us how external recruitment affects transitions 

between each state. The external recruitment matrix requires only one parameter a, 
being the probability that an empty patch is colonised by external recruits in a time 
step. Therefore we have that 

/(I - af 2a(l - a) a^X 
R = i 1 - a a . (3) 

V 1/ 

The new transition matrix A, with catastrophes, colonisation and external recruit- 
ment combined is given hy A = E x C x R. This order of matrix multiplication 
assumes that external recruitment occurs last, after the extinction and colonisation 
processes have occurred. The order of matrix multiplication does influence the results 
in this case, since we are interested in the first eigenvector (instead of the eigenvalue) 
of the matrix A. However, the order does not appear to have too much of an influence 
on the results from our investigations. This is partly due to the fact that we are only 
interested in finding the distance that minimises the probability of being extinct in 
equilibrium, rather than the actual probabilities of being in each of the states. The 
stationary distribution of A, which we obtain from the first eigenvector, tells us the 
probabilities of being in each of the states when the system is in equilibrium (Grim- 
mett and Stirzaker (1992)). From the stationary distribution we can find an analytic 
solution for the placement of marine reserves in a network, by examining the first ele- 
ment within the stationary distribution; by taking one minus the first element within 
the stationary distribution we can calculate the probability that a species is extant in 
equilibrium. This is the probability which we need to maximize to effectively use this 
model to manage marine reserves. 

The graph of one minus the first element in the first eigenvector is presented below 
(Figure (8)). We can see that when the recruitment probability is small, we find similar 
results to those for the baseline model. The optimal inter-reserve distance does not 
appear to change too much following the introduction of external recruitment, with 
only a slight increase found. The optimal intcr-rcscrvc distance also appears to be 
relatively independent of the external recruitment probability. 

5 Model with Local Extinction 

Another important aspect of metapopulations we consider in this investigation is the 
process of local extinction caused by non-catastrophic events (Akcakaya and Ginzburg 
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Figure 8: The probability of survival with external recruitment; r — 0.5, a — 0.1, /j, — 5 
and a = 0.05 (solid line), 0.075 (dots) 0.10 (dashed line), as three different local recruitment 
probabilities. 



(1991)). We construct a transition matrix L which contains the probabihties of tran- 
sitions between states due solely to local extinction events (e.g. which may occur 
through unfavourable local environmental conditions). The construction of L is sim- 
ilar to that for the matrix R; the matrix L only requires one parameter b, being the 
probability of local extinction within a time step. Therefore we have that 



The new transition matrix A, with catastrophes, colonisation, external recruitment 
and local extinction is given hy A = ExLxCxR. After constructing this Markov 
chain we perform the same analysis as with the previous model and find the spacing 
between reserves that maximizes the viability of the population (Figure (9)). 

It can be seen that local patch extinction does not have too much of an effect 
on optimal inter-reserve spacing. The probability of being extant in equilibrium does 
however decrease, and thus the determination of the optimal inter-reserve distance 
becomes of more significance for ensuring maximum population persistence. 
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Figure 9: The probability of survival with local extinction; r = 0.5, a. — 0.1, = 5, a = 0.1 
and h = 0.025 (sohd hne), 0.05 (dots) and 0.1 (dashed line), as three local extinction 
probabilities. 

6 Discussion 

Ecologists and marine park managers have for a great deal of time relied upon practical 
experience and qualitative methods to design marine reserves for conservation purposes 
(Shafer (2001)). The additional inclusion of a variety of stake holders into the process 
of marine reserve design may also detract from the formulation of a conservation 
management plan that ensures the survival of a particular species. We have provided 
methods that allow the precise calculation of the optimal inter-reserve distance. This 
allows decision makers to have quick estimates for optimal inter-reserve spacing, and 
thus quantitative information for discussion with stake holders. 

The results for our baseline model showed how to optimally space marine reserves 
under the competing processes of colonisation and catastrophe. This analysis involved 
maximising the second eigenvalue of the transition probability matrix, which results 
in the longest possible persistence of the species. We found that when the mean 
catastrophe size is smaller than, or comparable to, the mean dispersal distance, there 
exists a unique optimal inter-reserve spacing. As the mean catastrophe size increases, 
we arrive at a continuum of optimal inter-reserve distances. In this latter situation, 
the minimum of these optimal distances should usually be selected, to provide the 
maximum probability of colonisation. 

We then included the possibility of external recruitment and local extinction events. 
With the expansion of the model there is now no absorbing state and the first element 
of the first eigenvector tells us the equilibrium probability that the species is absent. 
To calculate the optimal reserve design strategy we minimize this probability. We 
found similar behaviour to the situation where there is no external recruitment and 
local extinction, but the importance of determining the optimal inter-reserve spacing 
was highlighted. 

As more information becomes available concerning species occupying marine land- 
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scapes, more accurate predictions for marine reserve design will be achievable. The 
values we have adopted in this paper appear typical for many species, with larval 
dispersal distances in the range of 10 — 100 km (Cowen et al. (2006)). 

Our results are of most interest in situations were the process of dispersal and 
catastrophe are on a similar scale. The extension of our results to the multiple-reserve 
network will also be interesting, as catastrophe and colonisation will always be on a 
similar scale, at least for some subset of the whole network. In such situations it will 
also be of interest to look at combining our results with those of others concerning 
marine reserve design, sec for example McCarthy ct al. (2005, 2006). In such situations 
there will always exist a trade-off between the size of the reserves, the number of 
reserves and the spacing, and the interaction between each with catastrophes will be 
of great importance (Airame et al. (2003)). 
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